
#ifndef __CRYSTAL_H_
#define __CRYSTAL_H_
#include <list>
#include <space.h>
#include <vector.h>
#include <exception.h>

using namespace std;

struct Atom{
public:
	unsigned int id;
	double radius;

	//relative coordinate of this atom in a lattice.
	COORDD dPos;
	Space* sp;

	Atom(Vector<double> p, double r, Space* sp1);

	//get absolute coordinate of this atom.
	__inline__ COORDD getAbsCoord(COORDD& ucCoord) const;
};


class UnitCell{
private:
	int get_num_atoms1()	;
public:

	//position of this unit cell.
	UNITCELLCOORD iPos;

	// number of atoms in this cell;
	unsigned int num_atoms;

	//the number of atoms, e.g. for a cubic lattice, a unit cell has "1" atom
	unsigned int num_atoms1;

	//all of the unit cell actually share the same list of atoms.
	//it points the atom list in the object of "Crystal"
	const ATOMLIST* atoms;
	const Space* sp;

	UnitCell(UNITCELLCOORD _iPos, const ATOMLIST* _atoms, const Space* _sp);

	COORDD getAbsCoord() const{
		return sp->getAbsCoord(iPos);
	}

	bool hasAtom(const ATOMCOORD& _co);

};


class Crystal{
public:

	const ATOMLIST* atoms;
	const UnitCell* cell;

	//space
	const Space* sp;

	Crystal(const Space* _sp, const ATOMLIST* _atoms, const UnitCell* _cell);

	__inline__ bool hasCell(const COORDI& cm) const;
	static bool hasCell(const COORDI& cm, const COORDI& dim);

	//get the neighbors of a cell.
	UNITCELLLIST getNeighbourCells(const COORDI& center);

	//get the neighbors coordinates of a cell.
	static UNITCELLCOORDLIST getNeighbourCoords(const COORDI& center, const DIMI& dim);

	//get the neighbors coordinates of a cell.
	UNITCELLCOORDLIST getNeighbourCoords(const COORDI& center);

	//build a lattice based on unit cells
	static Crystal* create();
};



/*************Member functions of struct Atom***********/
Atom::Atom(COORDD p, double r, Space* sp1){
	radius = r;
	sp = sp1;
	dPos = p;
}

__inline__ COORDD Atom::getAbsCoord(COORDD& ucCoord) const{
	return sp->getAbsCoordWithOffset(ucCoord, dPos);
}

/*************Member functions of struct UnitCell***********/


UnitCell::UnitCell(UNITCELLCOORD _iPos, const ATOMLIST* _atoms, const Space* _sp){
	iPos = _iPos;
	sp = _sp;
	atoms = _atoms;
	num_atoms = atoms->size();
	num_atoms1 = get_num_atoms1();
}

bool UnitCell::hasAtom(const ATOMCOORD& _co){

	ATOMLIST::const_iterator a;
	for(a=atoms->begin(); a!=atoms->end(); a++){
		ATOMCOORD ucCoord = sp->getAbsCoord(iPos);
		//if((_co - (sp->getAbsCoord(iPos) + a->dPos)).abs() < 0.001){ //if(abs.isSameWith(_co))
		//	return true;
		//}
		if(a->getAbsCoord(ucCoord).isSameWith(_co)) return true;
	}
	return false;
}

int UnitCell::get_num_atoms1(){

	if(atoms == NULL)
		THROW_LOGIC_EXCEPTION("Parameters much be set to UnitCell before calculating num_atoms1");

	bool hasXDirection = !(sp->getbasisX().isZero());
	bool hasYDirection = !(sp->getbasisY().isZero());
	bool hasZDirection = !(sp->getbasisZ().isZero());

	//check the number of dimensions
	int dim1 = (hasXDirection?1:0) + (hasYDirection?1:0) + (hasZDirection?1:0);

	//construct a lattice with dimension "crystalDim", and specify the center cell.
	DIMI crystalDim;
	COORDI centerCell;
	switch(dim1){
		case 1:
			if(hasXDirection){
			  crystalDim = DIMI(1, 1, 3);
			  centerCell = DIMI(0, 0, 1);
			}else if(hasYDirection){
			  crystalDim = DIMI(1, 3, 1);
			  centerCell = DIMI(0, 1, 0);
			}else if(hasZDirection){
			  crystalDim = DIMI(3, 1, 1);
			  centerCell = DIMI(1, 0, 0);
			}
			break;
		case 2:
			if(!hasXDirection){
				crystalDim = DIMI(3, 3, 1);
				centerCell = DIMI(1, 1, 0);
			}else if(!hasYDirection){
				crystalDim = DIMI(3, 1, 3);
				centerCell = DIMI(1, 0, 1);
			}else if(!hasZDirection){
				crystalDim = DIMI(1, 3, 3);
				centerCell = DIMI(0, 1, 1);
			}
			break;
		case 3:
			crystalDim = DIMI(3, 3, 3);
			centerCell = DIMI(1, 1, 1);
			break;
		default:

			break;
	}

	//central cell
	UnitCell cCell(*this);
	cCell.iPos = centerCell;

	//neighbor cell
	UnitCell nCell(*this);

	int* repCount = (int*)malloc(sizeof(int) * cCell.num_atoms);
	for(int i=0;i<cCell.num_atoms;i++)
		repCount[i]=1;

	UNITCELLCOORDLIST cs = Crystal::getNeighbourCoords(centerCell, crystalDim);
	UNITCELLCOORDLIST::const_iterator cii;

	ATOMLIST::const_iterator ai;

	COORDD ucCoord = sp->getAbsCoord(centerCell);

	cout<<"CenterCell:"<<"{"<<centerCell.x<<","<<centerCell.y<<","<<centerCell.z<<"}"<<endl;
	cout<<"CenterCell1:"<<"{"<<ucCoord.x<<","<<ucCoord.y<<","<<ucCoord.z<<"}"<<endl;

	for(cii=cs.begin(); cii != cs.end(); cii++){
		nCell.iPos = *cii;
		//cout<<"{"<<nCell.iPos.x<<","<<nCell.iPos.y<<","<<nCell.iPos.z<<"}"<<endl;
		int i=0;
		for(ai=cCell.atoms->begin(); ai != cCell.atoms->end(); ai++){
			if(nCell.hasAtom(ai->getAbsCoord(ucCoord))){
				repCount[i]++;
			}
			i++;
			//cout<<"{"<<ai->getAbsCoord(ucCoord).x<<","<<ai->getAbsCoord(ucCoord).y<<","<<ai->getAbsCoord(ucCoord).z<<"}"<<endl;
		}

		//cout<<endl;
	}

	//for(int i=0;i<cCell.num_atoms;i++)
	//	cout<<repCount[i]<<" "<<endl;

	long long temp = 1;
	long long temp1 = 0;
	for(unsigned int i=0;i<cCell.num_atoms;i++)
		temp *= repCount[i];
	for(unsigned int i=0;i<cCell.num_atoms;i++)
		temp1 += temp/repCount[i];

	cout<<"size:"<<cs.size()<<endl;
	cout<<"num_atoms1:"<<temp1/temp<<endl;

	free(repCount);
	return temp1/temp;
}

/*************Member functions of class Crystal***********/
Crystal::Crystal(const Space* _sp, const ATOMLIST* _atoms, const UnitCell* _cell){
	sp = _sp;
	atoms = _atoms;
	cell = _cell;
}

Crystal* Crystal::create(){
	//set up space
	Space* sp = new Space(DIMI(10, 1, 1), VECTORD(0.0, 0.0, 0.0),
			VECTORD(0.0, 0.0, 0.0), VECTORD(1.0, 0.0, 0.0));

	ATOMLIST* al = new ATOMLIST();

  //cubic
	al->push_back(Atom(COORDD(0.0, 0.0, 0.0), 0.5, sp));
	al->push_back(Atom(COORDD(0.0, 0.0, 1.0), 0.5, sp));
	al->push_back(Atom(COORDD(0.0, 1.0, 0.0), 0.5, sp));
	al->push_back(Atom(COORDD(0.0, 1.0, 1.0), 0.5, sp));
	al->push_back(Atom(COORDD(1.0, 0.0, 0.0), 0.5, sp));
	al->push_back(Atom(COORDD(1.0, 0.0, 1.0), 0.5, sp));
	al->push_back(Atom(COORDD(1.0, 1.0, 0.0), 0.5, sp));
	al->push_back(Atom(COORDD(1.0, 1.0, 1.0), 0.5, sp));

/*
	//bcc
	al->push_back(Atom(COORDD(0.0, 0.0, 0.0), 0.25, sp));
	al->push_back(Atom(COORDD(0.0, 0.0, 1.0), 0.25, sp));
	al->push_back(Atom(COORDD(0.0, 1.0, 0.0), 0.25, sp));
	al->push_back(Atom(COORDD(0.0, 1.0, 1.0), 0.25, sp));
	al->push_back(Atom(COORDD(1.0, 0.0, 0.0), 0.25, sp));
	al->push_back(Atom(COORDD(1.0, 0.0, 1.0), 0.25, sp));
	al->push_back(Atom(COORDD(1.0, 1.0, 0.0), 0.25, sp));
	al->push_back(Atom(COORDD(1.0, 1.0, 1.0), 0.25, sp));
	al->push_back(Atom(COORDD(0.5, 0.5, 0.5), 0.25, sp));
*/
	UnitCell* uc = new UnitCell(COORDI(0, 0, 0), al, sp);

	Crystal* c = new Crystal(sp, al, uc);
	return c;

}


//finish
bool Crystal::hasCell(const COORDI& co, const DIMI& dim) {
		 if(co.z<0 || co.z > dim.z-1 ||
			co.x<0 || co.x > dim.x-1 ||
			co.y<0 || co.y > dim.y-1)
		return false;
	return true;
}

//finish
bool Crystal::hasCell(const COORDI& co) const{
	return hasCell(co, this->sp->getDim());
}

//finish
UNITCELLCOORDLIST Crystal::getNeighbourCoords(const COORDI& center){
	return Crystal::getNeighbourCoords(center, this->sp->getDim());
}

//finish
UNITCELLCOORDLIST Crystal::getNeighbourCoords(const COORDI& center, const DIMI& dim){
	UNITCELLCOORDLIST ns;
	for(int zz=-1;zz<=1; zz++){
		for(int yy=-1; yy<=1; yy++){
			for(int xx=-1; xx<=1; xx++){
				UNITCELLCOORD newCell = UNITCELLCOORD(zz+center.z, yy+center.y, xx+center.x);
				//cout<<"{"<<newCell.x<<","<<newCell.y<<","<<newCell.z<<"}"<<endl;
		        //check if the new cell is in the crystal expand by dimension "dim"
				if(Crystal::hasCell(newCell, dim)
				&& (xx != 0 || yy!=0 || zz!=0)){
					ns.push_back(newCell);
					//cout<<"{"<<newCell.x<<","<<newCell.y<<","<<newCell.z<<"}"<<endl;
				}
			}
		}
	}
	return ns;
}

#endif
